#####################################################################
# Cook, Blas, Carroll, and Sinha 09/16/2016
# For "Two Wrongs Make a Right"  
# Results from MC experiments 
#####################################################################

setwd("C:/Users/sjcook/Downloads/dataverse_files") #Specify working directory 
rm(list = ls()) 

source("Experiment1.R") #Runs Experiment #1 (as described in the main text)
source("Experiment2.R") #Runs Experiment #2 (as described in the main text)

###keep only results for the coefficent estimates, standard error, and coverage for beta_0 and beta_1 from both experiments####
rm(list=setdiff(ls(), c("b0_est_ex1", "b1_est_ex1","se.b0_ex1","se.b1_ex1","cover.b0_ex1","cover.b1_ex1","b0_est_ex2", "b1_est_ex2","se.b0_ex2","se.b1_ex2","cover.b0_ex2","cover.b1_ex2")))

################################################
#Calculate Estimator Comparison Criteria 
################################################

###true values for beta_0 and beta_1###
beta0_true = -1 
beta1_true =  1

#################
#Experiment 1 
#################
bias.probit.b0_ex1 <-apply(beta0_true-b0_est_ex1,2,mean,na.rm=TRUE)
bias.probit.b1_ex1 <-apply(beta1_true-b1_est_ex1,2,mean,na.rm=TRUE)

st.dev.b0_ex1 <-apply(beta0_true-b0_est_ex1,2,sd,na.rm=TRUE)
st.dev.b1_ex1 <-apply(beta1_true-b1_est_ex1,2,sd,na.rm=TRUE)

se.probit.b0_ex1 <-apply(se.b0_ex1,2,mean,na.rm=TRUE)
se.probit.b1_ex1 <-apply(se.b1_ex1,2,mean,na.rm=TRUE)

mse.probit.b0_ex1 <-apply((beta0_true-b0_est_ex1)^2,2,mean,na.rm=TRUE)
mse.probit.b1_ex1 <-apply((beta1_true-b1_est_ex1)^2,2,mean,na.rm=TRUE)

cover.probit.b0_ex1 <-apply(cover.b0_ex1,2,mean,na.rm=TRUE)*100
cover.probit.b1_ex1 <-apply(cover.b1_ex1,2,mean,na.rm=TRUE)*100

#################
#Experiment 2
#################
bias.probit.b0_ex2 <-apply(beta0_true-b0_est_ex2,2,mean,na.rm=TRUE)
bias.probit.b1_ex2 <-apply(beta1_true-b1_est_ex2,2,mean,na.rm=TRUE)

st.dev.b0_ex2 <-apply(beta0_true-b0_est_ex2,2,sd,na.rm=TRUE)
st.dev.b1_ex2 <-apply(beta1_true-b1_est_ex2,2,sd,na.rm=TRUE)

se.probit.b0_ex2 <-apply(se.b0_ex2,2,mean,na.rm=TRUE)
se.probit.b1_ex2 <-apply(se.b1_ex2,2,mean,na.rm=TRUE)

mse.probit.b0_ex2 <-apply((beta0_true-b0_est_ex2)^2,2,mean,na.rm=TRUE)
mse.probit.b1_ex2 <-apply((beta1_true-b1_est_ex2)^2,2,mean,na.rm=TRUE)

cover.probit.b0_ex2 <-apply(cover.b0_ex2,2,mean,na.rm=TRUE)*100
cover.probit.b1_ex2 <-apply(cover.b1_ex2,2,mean,na.rm=TRUE)*100


################################################
# Caculate Marginal Effects 
################################################
#Experiment 1 
#################
me.true <-pnorm(beta0_true+beta1_true)-pnorm(beta0_true)

p11_ex1 <-apply(b0_est_ex1, 2,pnorm)
p12_ex1 <-apply(b0_est_ex1+b1_est_ex1, 2,pnorm)

me.est_ex1 <-p12_ex1-p11_ex1

bias.me_ex1 <-round(apply(me.est_ex1,2,function(x){mean(x-me.true,na.rm=TRUE)}),3)
se.me_ex1  <-round(sqrt(apply(me.est_ex1,2,var,na.rm=TRUE)),3)
mse.me_ex1 <-round(apply(me.est_ex1,2,function(x){mean((x-me.true)^2,na.rm=TRUE)}),3)

#################
#Experiment 2
#################

p11_ex2<-apply(b0_est_ex2, 2,pnorm)
p12_ex2<-apply(b0_est_ex2+b1_est_ex2, 2,pnorm)

me.est_ex2 <-p12_ex2-p11_ex2

bias.me_ex2 <-round(apply(me.est_ex2,2,function(x){mean(x-me.true,na.rm=TRUE)}),3)
se.me_ex2  <-round(sqrt(apply(me.est_ex2,2,var,na.rm=TRUE)),3)
mse.me_ex2 <-round(apply(me.est_ex2,2,function(x){mean((x-me.true)^2,na.rm=TRUE)}),3)





################################################
#Table 1: Simulation Study Results 
################################################
#Experiment 1 
#################
# Beta_0
round(rbind(bias.probit.b0_ex1,st.dev.b0_ex1,se.probit.b0_ex1,mse.probit.b0_ex1,cover.probit.b0_ex1),3)
# Beta_1
round(rbind(bias.probit.b1_ex1,st.dev.b1_ex1,se.probit.b1_ex1,mse.probit.b1_ex1,cover.probit.b1_ex1),3)
#################
#Experiment 2
#################
# Beta_0
round(rbind(bias.probit.b0_ex2,st.dev.b0_ex2,se.probit.b0_ex2,mse.probit.b0_ex2,cover.probit.b0_ex2),3)
# Beta_1
round(rbind(bias.probit.b1_ex2,st.dev.b1_ex2,se.probit.b1_ex2,mse.probit.b1_ex2,cover.probit.b1_ex2),3)


################################################
#Table 2: Marginal Effects 
################################################
#Experiment 1 
#################
round(rbind(bias.me_ex1,se.me_ex1,mse.me_ex1),3)  
#################
#Experiment 2
#################
round(rbind(bias.me_ex2,se.me_ex2,mse.me_ex2),3)  
